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ABSTRACT 

We present a general formalism for computing self-consistent, numerical solutions to 
the time-dependent radiative transfer equation in low velocity, multi-level ions under- 
going radiative interactions. Recent studies of time-dependent radiative transfer have 
focused on radiation hydrodynamic and magnetohydrodynamic effects without lines, 
or have solved time-independent equations for the radiation field simultaneously with 
time-dependent equations for the state of the medium. In this paper, we provide a fully 
time-dependent numerical solution to the radiative transfer and atomic rate equations 
for a medium irradiated by an external source of photons. We use Accelerated Lambda 
Iteration to achieve convergence of the radiation field and atomic states. We perform 
calculations for a three-level atomic model that illustrates important time-dependent 
effects. We demonstrate that our method provides an efficient, accurate solution to 
the time-dependent radiative transfer problem. Finally, we characterize astrophysical 
scenarios in which we expect our solutions to be important. 
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1 INTRODUCTION 

The treatment of radiation transfer in complex, dynamic 
physical systems is crucial for theoretical modelling in a 
wide variety of astrophysical contexts. In this paper, we ex- 
plore time-dependent effects in line transfer for low velocity 
media illuminated by external sources. Examples of such 
systems are found in absorbing material irradiated by light 
from gamma-ray bursts (GRBs) and active galactic nuclei 
(AGNs). In absorbing media, external photons enter at the 
boundary of the system, and are subsequently redistributed 
in angle and frequency by atomic interactions in the inte- 
rior. This leads to a complicated, time-dependent coupling 
of the radiation field and the state of the medium. 

Previous studies have investigated time-dependent ra- 
diative transfer in a number of regimes, under differ- 
ent approximations. Several papers have explored steady- 
state, non-LTE solutions to the radiative transfer equa- 
tion (RTE) for time-dependent media. These works were 
done in th e con te xts of magnetohy drodynamics (e.g., 
iHavek et al] l2010l: iDavis et all bOld) radiation hydro- 



dynamics (e.g., Krumholz et alj 20121). smoo thed parti- 
cle hydrodynamics 



fe.g.. iRundle et all [2OI0I). transport 



of cosmological ionization fronts (e.g., Whalen fc NormanI 



2006I1. molecu l ar ba nds in planetary atmospheres (e.( 
Kutepov et al. 199ll) , line transfer in stellar atmospheres 



(e.g., Sellmaier et al.l 19931). multi- level tran sfer in moving 



media (e.g.. Illauschildt ]l993l: iHarpedf 1994I ). atmospheres 
of hot stars (e.g., Hube nv fc Lanall995l ). clumpy molecular 



Email: miva3@mail.gatech.edu 



clouds (e.g.. iJuvela 1997i), and line e mission from interstellar 
clouds (e.g.. lJuvela k, Padoan|[2005[ ). 

A few authors have investigated solutions to the 
time-dependent radiative tr ansfer equa ti on (T DRTE). 
iPerna fc Lazzatil (|2002l ') and iPerna et abl l|200a ) treated 
time-dependent radiative transfer in static, dusty me- 
dia in the optically th i n lim it, wi thout radiation s ource 
terms. iGnedin fc Abell I2001I and lAbel fc Waiidel^ |2002| 
also solved the time-dependent radiative transfer equation 
for cosmological r eioniz ation, neglecting photon scattering. 
iMihalas fc Kleii] (|l982h developed a 'mixed-frame' formal- 
ism for solving the time-dependent equations in the context 
of fluid fl ow, and investigated numerical methods for their 
solution. Kunasj lfl98S ) calculated the numerical solution to 
the time-dependent transfer equation for resonance scatter- 
ing in a static medium consisting of two-level atoms. 

Recently, iHubenv fc Burrow^ (|2007h produced a self- 
consistent solution to the TDRTE for a time-dependent, 
non-LTE medium in radiation hydrodynamics. Their calcu- 
lation focused on the propagation of neutrinos in supernova 
explosions, and did not explicitly treat line transfer of pho- 
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tons (see also lAbdikamalov et al.|[2012l . for further discussion 
of time-dependent neutrino radiation transport). 

In this paper, we neglect fluid flow, and focus instead on 
time-dependent effects in the line transfer of photons. In ex- 
ternally driven absorbing material, the time derivative of the 
RTE cannot be neglected due to the variability of the incom- 
ing radiation at the boundary of the system. The resulting 
solution depends on several crucial timescales: the variabil- 
ity timescale of the external source, Tvr; the time interval 
between characteristic radiative interactions, T^of; and the 
light travel time across the width of the medium, Tlc- For 
optically thick systems, T^ei <S TLc- If Tvr <^ Tref , we argue 
that the source term for atomic interactions is unimportant 
for solving the TDRTE, and a simple, approximate solution 
is appropriate. Conversely, if Tvr ^ TLc, then the RTE can 
be solved in steady-state, while the medium is treated as 
time-dependent; this situation has been discussed exhaus- 
tively in the literature. Thus, we are primarily interested in 
the 'intermediate regime', Tvr > Tref, TVr < TLc, and will 
perform a detailed analysis of this scenario in the work that 
follows. 

The solution to the non-LTE radiative transfer prob- 
lem has an enormous associated literature, a proper re- 
view of which is beyond the scope of this paper. In the 
work that follows, we extend the Accelerated Lambda Iter- 
ation (ALI) method to the case of time-dependent radiative 
transfer and atomic rate equations. The ALI approach em- 
ploys a computationally inexpensive approximation to the 
RTE formal solution in the rate equations for the atomic 
level populations. This technique has its ori gins in the 'core 
saturation' method of iRvbickil 11972'. 'l984|), and Cannon s 
operat or splitting procedure (fCannorii ,1973h . lOlson et al.l 
l|l986l) showed that the diagonal part of the operator char- 
acterizing the formal solution to the RTE provides an ef- 
ficient approximate operator for implementing ALI. Subse- 
quently, iRvbicki fc Hummed (|l99ll . 119921 ) developed a pow- 
erful formalism for applying ALI to non-LTE problems with 
multi-level ato mic models. For a detailed review of ALI, see 
lHubenvll|2003t ). 

In this paper, we extend the TDRTE solution o f lKunasd 

to time-dependent media with multi-level at oms, 
making extensive use of the iRvbicki fc Hummed l| 19921 ) for- 
malism. To date, we have found no previous work in the lit- 
erature that treats both the RTE and atomic rate equations 
as time variable for multi-level line transfer. Our method 
applies to general atomic models, including bound-bound 
and bound-free radiative transitions, coUisional interactions, 
electron scattering, and background processes such as free- 
free absorption and emission. We present calculations for a 
specific three-level atomic model that illustrates important 
time-dependent radiative transfer effects. We find that these 
effects are important for systems in the intermediate regime 
and show that the use of approximate solutions that neglect 
them can lead to significant errors when interpreting spec- 
troscopic observations. We consider the implications of our 
results for astrophysical systems, and provide a set of criteria 
for determining when the use of our formalism is necessary. 

Our paper has the following outline: in ij2] we describe 
our models for the atomic physics and external radiation 
sources used in our calculations; in 331 we review the basic 
equations for the TDRTE and atomic rate equations; in ^ 
we develop a general formalism for implementing the ALI 



solution; in ij5l we show the results of our calculations; and 
in Sj6l we discuss implications of our results for astrophysical 
systems. In Appendix|Xl we extend our technique to include 
implicit Runge-Kutta time integration of our equations. 



2 MODELS FOR ATOMIC PHYSICS AND 
EXTERNAL RADIATION 

We consider a plane parallel, finite slab of material that is 
illuminated by an external source of unpolarized radiation. 
The properties of the medium vary along the z axis. The 
cosine of the polar angle, 6, measured from the z axis, is 
denoted by /i. The system is assumed to be symmetric in 
the azimuthal direction. 

The external radiation propagates toward the slab with 
6 > 7r/2, reaching the medium at an arbitrary coordinate, 
Zo- The width of the slab along the z direction is denoted 
by W; thus, radiation exits the slab at Zmin = Zo — W. 
After interacting with the medium, the radiation field in 
the slab consists of rays with — tt ^ 9 ^ tt. We assume that 
an observer is positioned in the 9 — tt direction. 

To illustrate key properties of the time-dependent so- 
lutions, we perform calculations for a specific three-level 
atomic model. The atomic levels are labelled by j — 0, 1, 2, 
and are connected by bound-bound radiative transitions be- 
tween j = ^ 1, j = 2, and j = 1 2. We neglect 
other types of radiative and coUisional interactions in our 
calculations, though they are straightforwardly included in 
our formalism, which is completely general. 

The number density of the atoms is denoted by no. We 
define a reference transition for photoexcitation between the 
atomic levels j = — >■ 2, at line centre, with all atoms in the 
ground state. If the Einstein coefficient for this transition 
is Bo2, and the line profile is due to Doppler broadening 
in complete redistribution, we obtain a reference extinction 
coefficient, 

_ hBo2no 

" 47r3/26/c' ^ ' 

where b = y/2kBTa/ma is the Doppler width for an atom 
with temperature Ta and mass nia (ignoring possible effects 
of microturbulence). We define a reference optical depth, 

dr = -Xrofdz, (2) 

such that t{zo) — 0, and the maximum optical depth is 
T(2min) = Tmax ~ XraiW . We uote that our definitions of z 
and r are opposite to those used in the stellar atmospheres 
literature, in which the optical depth increases with distance 
from the observer. 

If we define a reference time, 

Trc^ = (cxrof)"^ , (3) 

then the light crossing time for radiation passing through 
the slab is 

Tlc = W/c = r„,axr,.ef. (4) 

A natural unit of radiation intensity can be defined by 
equating the inverse of the photoexcitation rate for white, 
isotropic radiation exciting the reference transition, to the 
reference time: 

To = CXrcf/-B02. (5) 
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Table 1. Einstein coefficients for the transitions in our atomic 
model. A, -Bex, and -Bom are the coefficients for spontaneous emis- 
sion, photoexcitation, and stimulated emission, respectively. 



Transition 




Bcx/B02 


Bom/-B02 


0-^2 


1.4 X 10"' 


1 


1.25 


1^2 


4.9 X 10^ 


0.44 


0.44 


1 


0.18 


1.7 X 10-5 


2.1 X 10-5 



When the photoexcitation rate for the reference transition 
is equal to B02I0, the number density of line centre photons 
is equal to no. 

The Einstein coefficients for the transitions in our 
atomic model, using the units described above, are sum- 
marized in Table [T] The values of the coefficients have been 
chosen so that the atomic model contains a resonance line 
(j = 2), a strong transition between the excited states 
( j = 1 •<-> 2) , and a weak transition between the ground and 
first excited states {j — -H- 1). This arrangement insures 
that the second excited state (j = 2) remains relatively un- 
populated, while the transition j = ^ 2 — >■ 1 provides an 
efficient mechanism for populating the first excited stateQ 

In M to make comparisons to the previous work of 
lKunasjll983l ). we also consider a restricted version of the 
atomic model that includes only the single transition j = 
0^2. 

For the calculations in this paper, we consider several 
basic models for the external radiation arriving at r = 0. We 
assume that the angular variation of the incident radiation is 
either isotropic or highly coUimated at the r = boundary 
of the slab. In the latter case, photons arrive from the ex- 
ternal source in a narrow range of angles around the line of 
sight. In both cases, the spectrum of the external radiation 
is assumed to be white across all three atomic transitions. 

We employ a heuristic model for the time variation of 
the arriving photons. This model combines a transient in- 
crease in intensity with a constant background intensity: 



lT(t,yi) = |/b + 



j4cxp 



-ilog^ (tAo)/a^ 



/(m), (6) 



where 7^"' is the specific intensity of the incident radiation 
at r = 0, Jb is a constant background intensity, and A is the 
amplitude of the transient part of the incoming radiation. 
The variable part of equation ((6| is a Gaussian pulse in 
logt, where the parameters to and a determine the centring 
and width of the pulse, respectively. The angular variation 
/(/i) is set to unity for isotropic sources, and takes the value 
/(/i) = 1) for coUimated sources. Thus, if ^4 S> /b, the 

external source emits photons in a transient pulse, whereas if 
^ > 7b, the external radiation represents a persistent source 

with superimposed variability. 

T o make comparisons to the previous work of lKunas3 
j 19831 ). we also employ a simple 'ramp' model in which the 



^ The atomic data in Tabic [T] was chosen to mimic transitions 
from the ground state (^09/2) *nd first excited state (^£'7/2) to 
a common upper state (^^7/2) of the Fe II ion in fine structure 
splitting. Nevertheless, our three-level model provides a simple, 
general mechanism for populating an excited state by a resonant 
transition; this mechanism is widely applicable to many ions. 



external source increases linearly from zero at t/Trd = 
to /b over the time interval IQ'^Tr of, remaining constant 
thereafter. 



3 FUNDAMENTAL EQUATIONS 
3.1 Radiative Transfer Equation 

We restrict the range of angular values to ^ £ [0, 1], and 
write the specific intensity for the unpolarized radiation 
field as two distinct functions, for 9 ^ 7r/2 and 

I~{t, z, fi) for e > tt/2. With this notation, the TDRTE can 
be written: 



c dt dz 



(7) 



The extinction and emission coefficients can be cast in a 
general form for bound-bound and bound -free radiative in- 
teractions (c.f., iRvbicki fc Hummeilll99^ ). We denote the 
frequency- angle rate coefficient for spontaneous emission or 
radiative recombination as Uij , depending on whether I — >■ j 
represents a bound-bound or bound-free transition. Note 
that Uij = if / ^ j- Similarly, Vij denotes stimulated 
emission or recombination for I > j, and absorption or pho- 
toionization for I < j. For the familiar case of bound-bound 
transitions with complete redistribution in the line profile, 
the rate coefficients become: 



Kj 



hvij 



(8) 



(9) 



where vij is the line centre frequency of the transition, and 
ipij is the absorption line profile. Note that under the as- 
sumption of complete redistribution, tpij = ipji. Similar 
equations can be written for bound-free interactions, or 
for different approxi mations for the redistribution of the 
line p rofile (see, e.g.. iRvbicki fc Hummei1ll992l : lUitenbroekl 
l200lD . In this notation, the extinction and emission coeffi- 
cients become: 



where 



1,3 



Xij = VijUi - VjiUj 



(10) 

(11) 

(12) 



and Xc, ?7c represent background processes or electron scat- 
tering; these contributions are treated separately in the nu- 
merical method from the bound-bound and bound-free in- 
teractions (see @. The E symbols in equations (|10|) and (|ll|l 
denote double sums over the j and I indices for all transi- 
tions satisfying the indicated condition. The number density 
of atoms in level I is denoted by ni. 

Equation ^ must be supplemented by initial and 
boundary conditions for the specific intensity in each an- 
gular range. The boundary conditions take the form: 



/j" [t, Zmiu, At) = 0. 



(13) 
(14) 
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The initial conditions are determined by solving the time 
independent radiative transfer equation (TIRTE), which is 
obtained by setting dl^/dt — >■ in equation ([7]). Equa- 
tions (I13|l and (|14[) . evaluated at t/T^et = 0, are then used 
as boundary conditions for the TIRTE. The resulting solu- 
tion is used for the initial conditions, I^(t — 0,z,^j,). The 
TIRTE can be solved with an approach that is similar to 
the time-dependent method (see 

In developing numerical methods for solving th e RTE, it 
is use ful to define the Feautrier variables (see, e.g.. iMihalad 
1 19781 . and the references therein): 

= (/+ + /,;) /2, (15) 

V, = {1+ - I-) /2. (16) 

We define the optical depth along a ray path as: 

dTi, = -Xfdz/fi, (17) 

where t^(zo) = 0. A change in variables to u^, v^, and r^, 
results in the following form for the transfer equations: 



at V OTv 



dvi, 



(18) 
(19) 



where 5*1, = ^v/Xv is the source function. Similarly, the 
boundary conditions can be rewritten as: 

Uv{t,Zo,fJ.) = V^{t,Zo,ll) + /™'(t,/i), (20) 
U,y{t, Zmin, tJ-) = —Vi,{t, Zmin, m) • (21) 

Equations (|18p and (|19p , along with the boundary conditions 
(|20p and (|21|l . form the basis of our numerical solution to 
the RTE. 



3.2 Atomic Rate Equations 

The extinction and emission coefficients for the radiation 
field depend directly on the level populations of the various 
bound atomic states, as well as the degree of ionization of the 
medium. The rate equation for each level j can be written: 



(22) 



where Rij and Cij are the rate coefficients for radiative and 
coUisional transitions between levels / — >■ j, respectively. The 
sum E; is over all atomic levels I ^ j. Using the notation of 
i)3.1l we write the radiative rates as: 



Ri 



(23) 



where in the second equality we assume that the line profile 
for the transition is symmetric in fj,. The coUisional rate co- 
efficients are assumed to have a dependence on the electron 
number density, rie, and temperature, Te- For the purposes 
of our calculations, we assume that Ue and Te are known 
functions of t and z, or that they may be calculated sepa- 
rately from the atomic level populations, using a previous 
iteration in our numerical solution (see Q. 

Equation (|22|) must be supplemented by an initial 



condition for each level j. These are determined by solv- 
ing the time independent atomic rate equations, obtained 
from equation H22p by setting drij /dt — > 0, and using 
Uv{t = OjZjh) for the radiation field (see i|3.ip . Note 
that the steady-state calculations determining our initial 
conditions are equ i valent to the problem described by 
iRvbicki fc Hummer] (|l992l '). with an external source added 
to the boundary conditions. 

The primary difficulty in computing solutions for the 
radiation field and level populations are that equations H18p . 
(|19l) . and (|22p are coupled. Thus, the populations, Uj, must 
be determined simultaneously with the radiation field, u^. 
In <2] we describe a numerical method for calculating self 
consistent values for these variables. 



4 NUMERICAL METHOD 
4.1 Discretization Scheme 

We solve the RTE and atomic rate equations using numerical 
quadrature on discrete grids. For the reference optical depth, 
we define a set of D points, {t}^~q , that are equally spaced 
logarithmically over the interval [rmin, Tmax]- The number 
of spatial points D is chosen such that there is adequate 
resolution over the full range of optical depth for photons in 
the atomic transitions. In the models of ^ we set Tmin = 

10"* and Tniax = 10=^. 

For the time grid, we define a set of K points, {tk}^SQ ■ 
The appropriate spacing and interval for the time grid de- 
pends on the model used for the external radiation, as well as 
the values of the parameters in equation As discussed in 
JT] we are primarily interested in the intermediate regime for 
which TvR ^ Tref and Tvr ^ Tlc; we therefore choose val- 
ues of the parameters that reproduce this behaviour (see SJS}. 
When using equation © for the external radiation, we em- 
ploy a linearly spaced time grid spanning an interval equal 
to several light crossing times. When using the ramp model 
in our comparisons to previous work, we employ a logarith- 
mically spaced time grid that exceeds the light crossing time 
by more than an order of magnitude. The use of linear spac- 
ing for the time grid in the former case is important, as 
resolution at the smallest time variation scale is needed for 
each decade of time in the integration interval. 

For the frequency grid, we use the variable 



c = {iy-iyij)/Aiyij, 
Auij = uij {h/c) , 



(24) 
(25) 



where x is the number of Doppler widths, Avij, from line 
centre of the transition Z — >■ j. In iJSl we present results of 
calculations using Doppler and Voigt line profiles. For the 
Doppler case, we define a set of A'' points, {x^'I^Iq , spaced 
linearly in the interval [0,4]. For the Voigt case with pa- 
rameter a = 0.01, we augment the Doppler grid with a set 
of logarithmically spaced points out to a:: = 20. We assume 
that Ai/jj and a are the same for all transitions in our atomic 
model; we can therefore reuse the same frequency grid for 
each transition. It is straightforward to define a set of fre- 
quency integration weights, for numerical quadrature. 
In the calculations that follow, we use the extended trape- 
zoidal rule to define the integration weights, breaking the 
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integral into linear Doppler and logarithmic Voigt sections 
as necessary. 

It is standard to define a set of M angular grid points, 
{Pmim^o I such that each /x^ is an abscissa for a Gaussian 
quadratur e rule with co r responding int egration weight Wm 
(see, e.g., iMihalaj Il978l : [c anno nl ll985l ). It is often conve- 
nient, based on the form of the angular integrals, to choose 
the weights and abscissas for Gauss-Legendre quadrature. 
For an isotropic external source, this choice is adequate. 
However, for a highly collimated source, the angular vari- 
ation is such that the external radiation is non-zero for a 
narrow range of angles centred on 6^ = tt; this range is 
much smaller than the grid resolution of our models. In 
this case, we have found it useful to use Gauss-Lobatto 
quadrature, which includes the endpoint of the integration 
interval, p = 1, corresponding to radiation traveling in the 

= Tf direction for the radiation field J^. Recall that we 
restrict n £ [0, 1], and have defined /""'(t, /x) = I~(t,Zo,fi) 
to consist of external photons traveling in directions with 
6 > n/2. The value of the angle- averaged external inten- 
sity is fixed at J°'"(t) = i m)- We then set 
wol'^^it, ^o)/2 — Jl^^{t), where wq and are the inte- 
gration weight and abscissa associated with = 1 {9 = -n). 
These are the only non-zero values of 7°"' in the boundary 
condition, and can be used directly with equation ((20} . This 
procedure is equivalent to solving two sets of coupled equa- 
tions for the collimated and angularly redistributed compo- 
nents of the radiation field. 

It is convenient to use a single index that contains the 
combined frequency and angle dependences of the radiation 
variables. We define Q = N x M grid points with the formula 

1 = n + mN. The values of a function, g, on the grids defined 
above, are written gv„{tk,Td, Hm) = Okdi- We will use this 
notation extensively in the work that follows. 



4.2 Solution to the Radiative Transfer Equation 

To solve the pa rtial diff e rentia l equations (PDEs) in (|18|) and 
(fT9l) . we follow lKunas3 j 19831 ) and use a method of lines ap- 
proach, in which we replace the PDEs with ordinary differen- 
tial equations (ODEs) by discretizing the spatial derivatives. 
It will prove useful to augment our spatial grid by adding a 
set of D — 1 points, {Td+i/2}'d=o ^ ihaX are located between 
and equidistant from Td and Td+i - If equations (|18|) and (|19p 
are evaluated at and respectively, we obtain: 



dt 



cxdi — T Udi + Sdi] , (26) 



dVd+l/2,i _ rv _ 
- Jd+l/2,i — 



dt 



, Ud+l.i — Udi 1 /„-7\ 

cXd+i/2,i [ Ud+i/2,« j • (27) 



We use the following notation: 



At' 



Td+l,i — Tdi, 
Tdi — Td-l,i, 



A.Tdi = Td+l/2,i — 'Td-\l2,i 

= (Ar+ + Ar,l) 12. 



(28) 
(29) 

(30) 



Equations ([26} and ([27} represent a set of {2D~-4)Q coupled 
ODEs for d = 1, . . . , D — 2 and i = 0, . . . Q — 1. Equations 
for the variables at the boundaries will be considered below. 

The time derivatives can be integrated according to the 
'9 method', which parametrizes the solution in terms of the 
Backward Euler and Crank-Nicholson solutions: 

Ukdi = Uk-i,di + Aifc {9fkdz + fjfk-i^di) , (31) 

Vk,d+l/2,i = Vk-l,d+l/2,i + 



Atk {9fk,d+l/2,i + ^fk-l,d+l/2,i) 



(32) 



where Atk = tk ~ tk-i and fj = 1 ~ 6. The parameter 9 is 
defined such that equations (|3Hl and ([32} are integrated by 
the Backward Euler method for 9 = 1, and by the Crank- 
Nicholson method for 6 = 1/2. We explore more complicated 
Runge-Kutta integration schemes in Appendix 1X1 

Equation (|32p can be used to eliminate the variables 
'Vk,d±i/2,i in (|3ip . Substituting ([32} into ([31} and using the 
explicit forms for in ([26)) - ([27} yields the following sys- 
tem of equationsjj 



-AkdiUk,d-l,i + BkdiUkdi — CkdiUk,d+l,i 
= 9Skdi + 'fjkdiSk-l,di + FkdiUk-l,d-l,i 
+ GkdiUk-l,di + HkdiUk-l,d+l,i 
+ TkdiVk-l,d+l/2,i — Tf.diVk-l,d-l/2,i, 



where 



{Sk,d--L/2,i + 9) Ar^^.Azkdi' 



Ckdi 



(33) 



(34) 



(35) 



{Sk,d+i/2,t + 9) Ar+^^ATkdi ' 

Bkdi = Skdi + 9 + Akdi + Ckdi, (36) 



and 



Fkdi 



^'7fc,d-l/2,i 



{Sk.d-l/2,i + 9) ATi^_^^^.ATkdi 
(^Vk.d+l/2,i 



(5fc,d+i/2,i + 9) Ar+_^_^. Arfcdi ' 

Gkdi = Skdi — Vkdi ~ Fkdi ~ Gkdi, 

fjkdi , 9 5k^d±l/2,i — Vk,d±l/2,i 



kdi 



ATk- 



+ 



ATkdi Sk,d±l/2,i+> 



(37) 

(38) 
(39) 
(40) 
(41) 



using the definitions 

Skdi = 1/ (cxkdiAtk) , (42) 
Sk,d±i/2,i = 1/ {cXk,d±-^/2,iAtk) , (43) 

Vkdi = V {Xk-l,di/Xkdi) , (44) 
^,^±1/2,1 = V {Xk-l,d±l/2,i/Xk,d±l/2,i) ■ (45) 

Equation ([33} represents a set of D — 2 equations for 
each value of k and i. To complete the system of equa- 
tions, we spatially discretize ([19} to first order at d = 



Our expre ssions di f fer slig htly from the corresponding equations 
(9a)-(10) in lKunasj l ll983l'l . We can recover tlie equations in the 
earlier work by setting fjkdiTVk.d±l/2,i ~^ V ^rid substituting the 
explicit form for the two-level source function in our formulas. 
Numerical tests showed no significant difference in the solutions 
using either form for the two-level atom. 
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and d = D — 1, and use the boundary conditions (|20|) - (|2ip 
to eliminate the variables Vk,d=o,i s^nd Vk.d^D-i,i- If the re- 
sulting expressions are integrated with respect to time as 
described above, we obtain: 



{Sk,d^O,i + + Ck,d=0,i) Uk,d=0,i — Ck ,d=0,iUk ,d-l ,i 
~ {5k,d^0,i — flk,d=0,i ~ Hk,d^O,i) Ufc-l,d=0,i 
— Hk,d^0,i1l-k-l,d=l,i + {3k,d=0,i + ^ IkT 
+ (rik,d=0,i — Sk,d=0,i) Jfe-l,i, 



(46) 



and 



— Ak,d^D~l,i'U'k,D-2,i 
+ {Sk,d=D~l,i + 9+Ak,d=D-l,i)uk,d=D-l,i 
= Gk,d=D~l,iUk-l.D-2,i 
+ {5k,D-l,i — 'r]k,D-l,i—Gk,d=D~l,i)'U-k~l,D-l,i, 

where 

a 

C, 



A 



k,d=0,i 
Hk,d=0,i 
k,d^D-lA 



Ar, 



k,d^O,i 
f}k,d=0,i 



Art 



Ar, 



Gk 



k,d=D-l,i 



k,D-l,i 
rjk.D-l.i 



Ar,: 



(47) 

(48) 
(49) 
(50) 
(51) 

' fc-l,I5-l,i 

Thus, if the radiation field is known for all values of d and 
i at fc — 1, then equations (|33p . (|46p . and (I47p represent Q 
tridiagonal, linear systems for the spatial variation of the 
radiation variables, Ukdi, at time tk- In matrix form, the 
systems of equations can be written: 

Lfei ■ Uki — 6Ski + Sk~i,i + Bfci, (52) 

where Uki, Ski, and Sk-i,i are D x 1 column vectors with 
components Ukdi, Skdi, and f]kdiSk-i.di, respectively. The 
symbol liki denotes a tridiagonal matrix, while the column 
vector Bfci depends on the values of the variables at the 
previous time step, as well as the boundary conditions; the 
components of L^i and Uki are directly determined by equa- 
tions IMl), dMI, and (f47| . 

The initial condition for the radiation field, Uk=o,di, is 
determined by the prescription of H3.ll the time indepen- 
dent RTE is solved using the value of the source function at 
tfc=o. The equations governing the time i ndependent eq ua- 
tion are well k nown (see, e.g., Chapter 6 o f lMihalaslll97a . or 
Appendix A of lRvbicki fc Hummei|[T99ll ). The formal solu- 
tion to the TDRTE can then be obtained from equation 1)52^ 
and the initial condition by solving the tridiagonal system 
for successive k = 1, . . . K — 1. The results can be structured 
as a block matrix equation: 



Ui = Ai • Si + Bi, 



(53) 



where Ai is a, K x K block lower-triangular matrix, the el- 
ements of which are D x D matrices, and Ui, Si, Bi are 
K X 1 block vectors, the components of which are D x 1 
column vectors. The block vector Bi has a complicated de- 
pendence on the initial and boundary conditions for the so- 
lution. The structure of the formal solution in terms of the 
source function values, Skdi, is also complicated, with each 
pair of (fe, d) indices coupling to other spatial grid points at 
previous times. 



A key property of equation (|52|l is that, when Uk-i.di 
and Vk-i,d±i/2,i are known, it can be solved as a tridiago- 
nal matrix equation of dimension D, independently for each 
frequency- angle index i. Using the line ar algebraic methods 
descri bed in Appendices A and B of iRvbicki fc Hummeij 
(|l99ll ), the values of Ukdi (and hence the formal solution) 
can be computed in 0{KDQ) operations. In addition, the 
diagonal elements of L^/ , which will be used for the simulta- 
neous solution to the RTE and rate equations in t|4.4l can be 
calculated simultaneously at little additional computational 
cost. 



4.3 Solution to Atomic Rate Equations 

The equations governing the atomic level populations in ((22} 
form a system of ODEs. Though the system exhibits no 
explicit coupling of the populations at different spatial grid 
points, coupling is introduced implicitly by the presence of 
the radiation field in the rate coefficients, Rij . 

The time derivative in equation (|22p can be integrated 
with the same technique used for the RTEs in i]4.2l Applying 
the 9 method, we obtain: 

nj^kd = nj^k-\,d + Atfe {Ofj,kd + Vfj^k-i,d) , (54) 

fj,kd = |^(C';j,ftd 4- Rlj,kd) ni^kd 



{Cjl^kd + Rjl,kd) '^i.fcdj • 



The radiative rate coefficients take the form: 
-1 



— ~7~ ~ [Ulj,i + Vlj^i Ukdi] ■ 



lj,kd ' 



(55) 



(56) 



where the integration weights for the double integral in equa- 
tion (|23|) are given by qi = hnWm- From equation H53p . the 
expression for Ukdi contains terms of the form [Ai ■ Si]^^, 
which couple various Skdi for different values of k and d. 
Because Skdi has a non-linear dependence on the level pop- 
ulations rij, equations (|54|l and (|55|l represent a non-linear 
system of equations for the populations, coupled in the in- 
dices j, k, and d. 

The initial conditions for the rij are obtained from equa- 
tion ((22} by setting the derivative on the left-hand side equal 
to zero, and using the initial condition for the radiation field, 
Uk=o,di in the rate coefficients (see H4.2|l . This results in a 
coupled non-linear system in the indices j and d. 

4.4 Simultaneous Solution with Accelerated 
Lambda Iteration 

The simplest simultaneous solution method for the RTE and 
atomic rate equations, often referred to as 'Lambda Itera- 
tion', can be outlined as follows: 

(i) Start with a set of values for the level populations at 
all grid points, denoted j.^. These populations determine 

the values Xkdi ^^'^ vldi equations (|lUp and (|lip . 

(ii) Using these values, compute the formal solution for 
the radiation field, ul^^ (where the symbol f indicates that 
the populations j,^ have been used to obtain the field). 

(iii) Substitute the field u\^^ into the linear systems of 
equation ((54} and solve for an updated set of populations 
rij^kd- 
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(iv) Check for convergence of the populations by forming 
the quantity: 



e = Max 



(57) 



where eA and en are the absolute and relative local error 
tolerances, respectively (chosen for the particular problem 
to be solved). The symbol Max indicates that the largest 
value for the quantity in parentheses, considered for each 
combination of j, k, and d, should be assigned to Q. 

(v) If 6 ^ 1, we consider the system of equations to have 
converged, representing a self consistent solution for the field 
and level populations. If B > 1, we replace nj^kd ' '^^ 
and repeat steps 1—4. 



j.kd^ 



Unfortunately, this scheme suffers from several deficien- 
cies, which render it unusable for many practical problems 
with Tmax S> 1. For detaile d discu s sions of these issues, see, 
for example. [Mihalad (|l978l ), lAueJ (|l99ll ). and the references 
therein. 

We can improve on the Lambda iteration scheme by 
using the ALI method. This technique replaces the formal 
solution of equation (|53|) with an approximate expression to 
be used in the iteration scheme: 



(58) 



where Aj is an appropriately chosen approximation to the 
full matrix, and the symbol f indicates that S|, Bt, and 
u| have been calculated using the quantities j.^. As the 
level populations converge, j,^ — > Uj^kd, and equation (|58p 
becomes identical to equation (153 



As originally shown by , Olson et al.l l| 19861 ). an efficient 
choice for A^ is the diagonal of the full matrix A^. With this 
choice, we can write equation (|58|l as: 



Ukdi = ^Afc 



Skdi 



(59) 



where AJ_ 



vj,^i = [J-'kii jjjj denotes the diagonal elements of the 
matrix Ai, which are calculated as described in M.2\ 

In principle, more complicated choices for A can lead 
to faster convergence. For example, instead of using the diag- 
onal of the full matrix, we could use the tridiagonal subma- 
trix of A . This leads to a more complicated expression than 
H59p. involving elements that cannot be obtained easily from 
^ki' thus, considerable computational effort is required to 
use the ALI method in this case. This situation is well known 
for multidimensional radiative transfer problems (see, e.g., 
iKunasz fc Qlsonlll988l : lAuer et al.|[l99i ). 

To implement the ALI method, we follow 
iRvbicki fc Hummeij |l992) and use an alternate for- 
mulation of equation (|59l) : 



Ukdi = O^tai (vkdi - 7?Li) + "L 



(60) 



can be approximated using the results from a previous it- 
eration. Therefore, the quanti ties r?e cancel in equation ((60} 
(c.f., the discussion in §2.1 of lRvbicki fc Hummer! [l992i ). 
The new matrix '^kdi is related to A^^j^ through 



"^kdi — Afcdi/xL 



(61) 



The f symbol appears in this equation because the approx- 
imate matrices are themselves constructed using the quan- 
tities nt j,^. The advantage of using '^kdi is that Ukdi now 
depends linearly on the updated quantities Uj^kd through 
r]kdi, whereas Skdi exhibits a non-linear dependence due to 
the factor of Xkdi in the denominator. 

The combination of equations (|56|l . and (|60p yields: 



Rlj,kdni,kd — Rjl,kdnj^kd 

= -r- — S Ulj^iUl^kd — Ujl^iUj^kd + Xlj,kdi 

i^O I 



l',3' 



(62) 



Non-linearities occur in equat ion (1621) due to terms of 
the f orm xij,kdinii ^kd- Following iRvbicki fc Hummeij (|l99ll . 
Il992l ). we linearize the equation by making the replacement 
Xi3,kdinii ^kd x\j,kdi^i' .kd^ This procedure results in the 
'preconditioned' expressions: 

Rlj,kd'ni^kd — Rjl,kd'nj,kd 
-1 



4-K 



E 

i=0 



UljAlll^kd — Ujl^iUj^kd + d^kdi 



(63) 



X E Ui'j'.i{x\j,kdini'M - Xij,kdin\, ,^^) 
I'.i' 



+ Xlj,kdiu\.^^ 



We can rewrite equation (|63|) using our RTE formal solution 
and rate equation discretization. Employing identity (|61|l 
and collecting terms about rii^kd, nj^kd, and riii ^kd, we obtain: 



where 



(64) 



fj,kd = Y2 ^{^Ij.kd + 'R'lj,kd)ni^kd 
I 

— {Cjl^kd + 'Tl.jl,kd)'nj^kd 



Tll,,kd = E ^ {^'^-.^ [l - {x]ukd^lxld}) 'Q^ld\ 



(65) 



+ 



and 



7^ 



lj,kd 



47r 

T 



Vlj,' ['^ld^ - S^*kd^ (S'L, - Vl,kd^/Xld^)] } 

Q-1 

E 77- (xlkdJxL) OKd^ E '^^'r.^- (66) 



kd ; ' "-kdi- 



Unlike the line transitions, the background and scattering 
contributions to r] are treated as quantities with values that 
are either fixed in the problem under consideration, or that 



3 As noted by IRvbicki fc Hummeil l ll99ll) . the alternative proce- 
dure Xljykdi'f^i' ,kd ~^ Xlj,kdi^l/ kd ^ecovers the Lambda Iteration 
scheme. 
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Substitution of equations (|64p - (|66p into ((54} yields explicit 
expressions for the linear systems. 

Lambda Iteration can now be replaced by an improved 
ALI scheme that exhibits faster convergence. The algorithm 
is essentially the same as that for Lambda Iteration except 
that the modified linear system is used to compute the up- 
dated level populations. 

The ALI rate equations have some attractive proper- 
ties. Given the initial conditions nj^k^o,d, Uk=o,di, and tak- 
ing Afc^Q^j — 0, the ALI iteration can be applied to succes- 
sive k=l,...K-l, using the values A^.^^^^, xLi.di, and 
^k-i di from the previous time step to initialize the ALI iter- 
ation in the next time step. If L is the number of levels in the 
atomic model, the solution of the rate equations at a given 
time step requires 0{DL^) operations, and each iteration 
requires 0{DQL'^) operations. Therefore, in the numerical 
implementation of our solution, the time steps form an outer 
loop, while the iteration over the level populations forms an 
inner loop. This design for our iterative solution is analogous 
to that used bv iHubeny fc B urrows (2007), who noted that 
employing level populations from the previous step to seed 
the iteration of the current step significantly increases the 
convergence rate in a time dependent calculation. Because 
the structure of our solution method is the same as in their 
work, our calculations can be straightforwardly implemented 
in applications requiring radiation hydrodynamics. 

It should be noted that one unattractive property of 
the ALI equations is that they introduce coupling between 
atomic levels that are not associated by transitions in the 
original rate equations, due to the terms containing R\j^kd- 
While equation (|64|) can be implemented directly with (|54|) . 
for many problems of interest substanti al simplifications can 
be achieved. As discussed in §§2.4-2.5 of lRvbicki fc Hummeij 
for atomic systems containing transitions that don't 
exhibit 'significant' frequency overlap, terms of the form 
Uvji^i can be neglected when Ij and I' j' denote dif- 
ferent transitions, by setting Rij^^d — >■ in equation (I64|) . 
The resulting linear system only couples levels that are as- 
sociated by transitions in the original rate equations. 

In the following, we exclusively consider systems that 
exhibit negligible frequency overlap in their transitions. For 
fixed Ij, only a single transition contributes to the extinction 
and emission, which are non-zero for a subset of the total 
frequency grid. In this range, xlj^kdi ~^ xLi. and the terms 
in (|65p becom^fl 

Tlij,kd = -r- — [f^ij,*^! ~ S^kdi) 

l^o""^ ^ ' (67) 

To obtain equation (I67p . we also neglected background and 
scattering processes in the extinction and emissivity. We 



* We note that the result I I67II can be obtained directly by sub- 
stituting equation I I58I I into the RHS of I I56I I for the case of 
non-overlapping lines in complete redistribution. No precondi- 
tioning of the equations is necessary due to a serendipitous can- 
cellation of the denominator in the source function. This does 
not hold, however, for mo re general cases (c.f., the discussion in 
iRvbicki fc HummejlToQlh . 



Table 2. Descriptions of the calculations presented in j5]of the 
text. The canonical model uses the three-level atom described 
in ^2] and the conditions for the external radiation set by equa- 
tion (|6)l. The parameter values for the latter are given in il5.2l 
The descriptions of the external source and atomic models in the 
third column denote deviations of the calculation from the canon- 
ical case. 



Model 


External Source 


Atomic Model 


Model I 


Linear increase to /b 


Two-level 


Model II 


Linear increase to 7b 


Canonical 


Model III 


Canonical 


Canonical 


Model IV 


Canonical 


j42i multiplied by ten 


Model V 


A//o = 1 


Canonical 


Model VI 


A//o = 1, /b//o = 0.1 


Canonical 



used this form of the ALI system to calculate the results 
presented in ij5] 



5 RESULTS 

We present results from several calculations using the atomic 
physics and external radiation models described in ij2] We 
performed a total of six calculations: Models I and II test 
our numerical solution against previous work in the litera- 
ture and explore the approach of our solutions toward an in- 
dependently computed steady-state solution. Model III rep- 
resents a canonical solution that exhibits a number of im- 
portant time- dependent effects associated with the transfer 
of radiation through the time-dependent medium. Finally, 
Models IV- VI explore the consequences of varying certain 
parameters that characterize the external radiation source 
and atomic physics. The essential model properties are sum- 
marized in Table [2j and are described in detail below. 

5.1 Comparison to Previous Work and Numerical 
Tests (Models I and II) 

As a ba s ic tes t of our code, we reproduced the results of 
iKunasd (119831 ) . using the methods of (Model I). The 
earlier work used a constant source for the external radi- 
ation and considered time-dependent radiative transfer in a 
static medium consisting of two-level atoms that effectively 
remained in the ground state. The line profile for the transi- 
tion was assumed to be a Doppler or Voigt profile; the latter 
used parameter a = 0.01. For a resonant transition in a two 
level system, a steady-state medium is a good approxima- 
tion, as each excitation event is followed by a spontaneous 
de-excitation to the ground state. We reproduced this effect 
with our time-dependent methods by restricting our atomic 
model to the j = — 2 transition. In this case, the 712 pop- 
ulation remains negligible compared to rio, but the source 
function for the transition is finite Q To compare with pre- 
vious work, we did not use the initial conditions described 
in ij ^4.2ll4^ but rather set all atoms in the ground state at 

t/Trcf = 0. 

We performed calculations using this restricted model 
for grid resolutions oi K ^ 170, D = 140, iV = 40, and 

^ That is, A2on-2/(7?02'^o — ^20^2)^ remains significant, though 
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M — 10. The reference depth points were spaced with 20 
points per decade in r. The temporal resolution, K, was 
chosen to approximately equal the resolution in the spatial 
grid. We used the ramp model for white, isotropic external 
radiation at the boundary. Our time grid was logarithmi- 
cally spaced over the interval [10"*, 2 x 10*] in units of T^cf- 
For the line profile, we used a Voigt function with parameter 
a = 0.01, and divided the frequency points evenly between 
the linear and logar i thmic sections of the grid (see tj4.1|l . In 
the work of lKunas 3 (|l983l ). the radiation field was separated 
into unscattered and diffuse components, consisting of pho- 
tons that underwent zero and at least one radiative interac- 
tion, respectively. We computed the unscattered component, 
junsc fixing the level populations at their converged values 
(which determined the extinction coefficient), and then set 
the source function to zero. The diffuse field was calculated 
as 7^'« = h- /r""- 

Figure [T] shows the results of our calculations, which 
can be compare d with the righ thand top and bottom pan- 
els of Fig. 1 in iKunasj l|l983h . The top panel shows the 
two-level source function for the diffuse component, J = 
(47r)~^ J^^ dQ J duip^I^^^, as a function of reference depth, 
T, at several times. The bottom panel shows the angle- 
averaged intensity, = {4-k)~^ J^^ dQI^^^ , at r = 0, as 
a function of x, for several times. The curves are all normal- 
ized to the maximum value of an independently calculated 
steady-state solution (see ij i|4.2fl473)) . The different curves 
correspond to monotonically increasing values of the time 
from bottom to top in each panel. These times repre sent the 
closest values on our grid to those listed in Table 1 o f lKunasd 
l|l983l l. To the extent we were able to compare, our results 
are in excellent agreement with the earher work, reproducing 
the main features and trends, and showing only slight quan- 
titative disagreement in one of the plots □ Though we show 
a comparison to only two figures for b revity, w e are a ble to 
reproduce all of the essential results of lKunas3 (| 19831 ). 

The main features of Figure [T] are propagation of the 
external radiation through the medium with increasing time 
in the top panel, and the approach of the curves toward 
an independently calculated steady state solution in both 
panels as i/Tref increases beyond the light crossing time Tlc- 

Using the same model for the external source, we per- 
formed a similar calculation using the three-level atomic 
model described in ^ and a Voigt line profile (Model II). 
In this case, the medium is no longer approximately static, 
leading to a non-negligible population in the first excited 
state, j ~ 1 (see tj5.2|l . This leads to a difference in the 
source functions for the j = — >■ 2 and j — 1 ~^ 2 tran- 
sitions compared to the two-level case, altering the diffuse 
radiation field in the medium. 

Figure [2] is identical to Figure [T] except that it shows 
results for the three-level model, using the radiation field 
within a; ^ 20 of line centre for the j — 2 transi- 
tion. With this atomic model, the radiation field exhibits 
a stronger depression at large reference depths, as well as 
a slightly softer line profile at r = 0. This is caused by the 
fact that, in contrast to the static two-level model, not every 



° In tlie righthand bottom panel of Figure 1 in iKunasj {Tgll), 
we obtain quantitative agreement if the labelling of curves 'b' — 
'k' is shifted down by one letter to 'a' - 'j'. 



Ill" 



111-' 



ur' 



iir» 




ur- 10"' 10° 10' 10- 111'' 




X 



Figure 1. Plot of the static, two- level source function, J for the 
J = — > 2 transition as a function of t (top panel) ; and the an- 
gle averaged intensity as a function of frequency in Doppler 
widths, X (bottom panel). In the bottom plot, the results are 
shown at T = 0. Calculations were performed for Model I (see 
i|5.1ll . For both plots, each curve represents a distinct time, in- 
creasing monotonically from bottom to top. The curves are all 
normalized to the maximum value of an independently calculated 
steady-state solution (see ii ^4.24j43)l . The figure may be c ompared 
with t he righthand top and bottom panels of Fig. 1 in iKunasa 
l|l983! ). 



excitation to j = 2 results in a decay to j = 0; a fraction of 
de-excitations are to the j — 1 level, determined by the val- 
ues of the populations and the Einstein coefficients for each 
transition. An important feature of Figures [1] and [2] is that 
the magnitude of the angle-averaged diffuse field becomes a 
significant fraction of the external Ji, only after an interval 
> Trof . This is expected, as Trcf represents the characteristic 
time interval between successive radiative interactions. This 
point is important for judging whether a given astrophysical 
system will exhibit time coupled effects in radiative transfer 
(see p. 

Figure [3] is identical to Figure [H except that it shows 
results at frequencies x ^ 20, measured from line centre of 
the j — 1 ^ 2 transition. The results are similar to those 
described above, except that the radiation experiences rela- 
tively smaller extinction at large reference depth. This is due 
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X 



Figure 2. Same as Figure [T] except for tlie tiiree level atomic 
model (Model II of i|5.1ll . The quantities shown were computed 
using the radiation field within x ^ 20 of line centre for the j = 
0—^2 transition. 



to the fact that the effective optical depth for the excited 
state transition is much smaller than that for the ground 
state resonant transition (see tj5.2p . As in Model I, the re- 
sults for Model II converge to an independently calculated 
steady-state as t bec omes larger th an Tlc- 

§IIId and §IV of lKunas3 l| 19831 ) present a detailed anal- 
ysis of the method of lines solution to the radiative transfer 
equation in a two-level, static medium, analysing both its 
stability and effectiveness for various test problems. We have 
come to similar conclusions regarding our solution for the 
three-level, time-dependent medium. Briefly, we note that 
for 9 < 0.6, unphysical oscillations can be introduced into 
the solution (though this can be mitigated somewhat by re- 
ducing the time and space grid spacings). Also, as the time 
and space grid resolution is increased, the results approach 
a final, converged solution with respect to the grid spacing. 
As in previous works, we find that distinct spacing in the 
time and space grids leads to an approximate representation 
of the propagation of the leading front of the external radi- 
ation. As long as the time interval for integration is not too 
large compared to TmaxTl-ef, we can use a linear time grid of 
sufficiently small spacing to provide an adequate sampling 
of the variation of the radiation field at all reference depths. 




iir' i()-' 111" 1(1' 1(1- Id'' 




X 



Figure 3. Same as Figure |2] except results are shown for the 
radiation field within x 20 of line centre for the j = 1 — > 2 
transition. 
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Figure 4. Maximum relative change in the level populations as 
a function of iteration number for the ALI and Lambda iteration 
methods applied to Model II of ^5.11 
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Apart from the solution to the RTE itself, the central 
concern of our method is the convergence of the atomic level 
populations. Figure [4] shows the maximum relative change 
in the populations of levels no, ni, and n2, for all time and 
space grid points, at each iteration during the solution. The 
curves marked by diamond and circle shaped symbols show 
the relative population change when Lambda and ALI it- 
eration is used, respectively. As expected, the ALI method 
has a steeper slope and converges more quickly for a given 
number of iterations than the Lambda method. The results 
shown in Figure [4] are fairly typical for the models presented 
in this paper; in fact, the convergence properties are often 
better for Models III-VI than for Model II. For practical im- 
plementation of our methods with more complicated atomic 
models, it will likely prove useful to supplement the ALI it- 
eration with additional acceleration methods, for example, 
Ng's method or the generalized minimum residual method 
(for a detai led discussion of acceleration techniques, see, e.g., 
lAueij|l99ll ). 

As we increased the spatial grid resolution, the number 
of required iterations to reach a given relative accuracy also 
increased; this is a well known phenomenon (see Olson e t al.l 
1 19861 ). However, this was typically not true of the time grid 
resolution. As discussed in tj4.4l the initial guess for the solu- 
tion at each time step is the converged solution for the pre- 
vious time step. The finer the time grid resolution, the more 
similar the solution is in consecutive steps. Thus, increasing 
the time grid resolution tends to reduce the number of iter- 
ations per step. Since a reduction in spatial grid spacing is 
typically accompanied by a reduction in time grid spacing, 
the two effects will tend to counteract one another, and the 
increase in spatial grid resolution will not adversely affect 
the required number of iterations as much as expected. This 
is an unexpected benefit from applying the ALI method to 
the time-dependent problem. 

5.2 Canonical Transient Pulse (Model III) 

A class of external sources of immediate interest for this 
paper are those exhibiting transient behaviour, in which the 
amplitude of the transient pulse, A, is much greater than the 
background intensity 7b. We constructed a canonical model 
for transient sources, using equation with parameter 
values Ib/Io = 10"^ A/Ia = 0.1, to/Trof = 5 x 10^ and 
a = 0.5 (Model III). We assumed that the external source 
was collimated, according to the prescription of tj4.1l In our 
calculations, we used grid resolutions of -R" = 280, D = 280, 
A'' — 20, and M — 40. The time grid was spaced linearly 
over the interval [0, 4rLc]- 

For our atomic model, we used the three-level system 
described in iJ21 with a common, Doppler line profile for all 
three transitions. The initial condition for the radiation field 
and atomic level populations was the steady-state solution 
at t/T,^i = (see ^ H4.2H4.3|) . 

We define the decrement in the radiation field associ- 
ated with a given transition as the value of u^, at line centre 
(x = 0) divided by the value in the line wing (x = 4). 
The decrement is denoted Uccnt /wwing . Figure [5] shows the 
results of our calculation for the the j — ^ 2 transi- 
tion, with direction /i = 1. The top panel shows the line 
decrement as a function of t at several reference depths, 
r = 0, 1,2. The time t is measured in units of Tiof. Larger 




Figure 5. Radiation field for the j = — > 2 transition with 
fi = I. The top panel shows the line decrement ticcnt /"wing (see 
the definition in i|5.2l l for several values of the reference depth. 
The bottom panel shows the radiation variable function of 

reference depth at several times. In the bottom panel, the variable 
Ui, is normalized to its maximum value in the t/T-^^{ = 1000 curve. 
The time is measured in units of T^^f. All calculations have been 
performed for the canonical Model III of ^5.21 



depths are not plotted, as the radiation field is negligible 
due to large extinction at r S> 1. The peak intensity occurs 
at t/Trcf ~ 500, when the external radiation has its maxi- 
mum value. At r = 0, there is little radiation re-emitted in 
the > direction, implying that « 7^/2 « 11"^^ / 2. The 
bottom panel of Figure [5] shows as a function of r for 
several times: i/T^cf = 1000, 1500,2000. These times corre- 
spond to Tlc, 1.5Tlc, and 2Tlc, respectively. The curves in 
the bottom panel are normalized to the maximum value of 
at t/Trcf = 1000. 
Figure [6] is similar to Figure (5] except that it shows 
results for the j = 1 ^ 2 transition. The radiation field is 
non-zero at larger reference depths than in the transition j — 
— > 2 because the effective optical depth for the j = 1 ^ 2 
transition is smaller. This is due to the reduced population 
of the excited state relative to the ground state, as well as 
the smaller Einstein absorption coefficient for the j = 1 — >■ 2 
transition compared to the resonant transition. It should be 
noted that most of the extinction in this transition occurs for 




Figure 6. Same as Figure O except results are shown for tiie 
j = 1 2 transition at reference optical depths r = 0, 100, 1000. 

T < 100. The radiation field at frequencies near the resonant 
transition, which is necessary for the excitation event j — >■ 
2 — ^ 1, is largely damped at r = 100. Thus, there are few 
excitation events to j = 1 at r > 100, and therefore little 
difference in the extinction from the j = 1 — !> 2 transition at 
r = 100 versus r = 1000. 

The smaller extinction is also evident in the bottom 
panel of Figure [SI in which the crossing of the spatial radi- 
ation profile through the r = Tmax boundary is evident as 
time increases. As expected, the maximum of the radiation 
pulse is shifted to later times at larger reference depths. This 
is due to the fact that the light crossing time is of order the 
variability timescale for the radiation source. This feature is 
absent from Figure \S\ because the reference depths plotted 
are small compared to TLc/Trof. 

The results for the radiation field at frequencies around 
line centre of the j = ^ 1 transition are similar to those 
shown in Figure [G] except for the absence of significant ex- 
tinction, due to the small value of the Einstein absorption 
coefficient for this transition. 

Figure [7] shows the angular profile of the radiation at 
several frequencies around line centre of the j = — 2 tran- 
sition. Results are plotted at a; = 0, 1, and 1.5. The top panel 
shows the values of the radiation field variable at r = 
for t/Trcf = SOOTrof, when the external source peaks at the 



Figure 7. Radiation variable 111/ ^ 

function of /i for several 
values of frequency in units of Doppler width, x. The top and 
bottom panels show results at r = and t = rmax for times 
t/Trcf = 500 and t/T,.^{ = 500 -|- Tmax, respectively. All calcula- 
tions were performed for the canonical Model III of ^5.2\ 



inner boundary. The bottom panel shows the values of Ui, 
at T — Tmax for t/Ticf ~ 1500 Trcf, wheu the peak of the 
external radiation passes through the outer boundary. We 
show the results for the diffuse field only; therefore, values 
at ^ = 1 are excluded from the plots. All curves are normal- 
ized to the maximum value of u at each frequency x. The 
plots demonstrate that the radiation field remains highly 
coUimated at the r = boundary but is largely isotropic 
at Tmax except at j: > 1. This result makes sense: because 
the optical depth in the j = — >■ 2 transition is large, pho- 
tons that reach r — Tmax have experienced many atomic 
interactions and have been redistributed isotropically in an- 
gle. By contrast, the angular distribution of photons at the 
T = boundary is largely set by the external radiation; 
only backscattered photons from nearby layers of the slab 
are redistributed in angle at t = 0. The medium therefore 
transitions from highly coUimated to largely isotropic with 
intermediate results at t > 1. 

Figure [8] is the same as Figure [7] except that it shows 
results for frequencies around the j — 1 ^ 2 transition. In 
both cases, the radiation with 6^ > is a small fraction of 
tt^ at T = 0. However, for the j = 1 — >■ 2 transition, the 
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Figure 8. Same as Figure [7] except for frequencies around the 
j = 1 2 transition. 



radiation field remains highly anisotropic even at t — Tmax- 
This is due to the fact that this transition has a small effec- 
tive optical depth (see below). Therefore, a small fraction of 
photons reaching r = Tmax have been absorbed and redis- 
tributed into angles fJ. ^ 1. The radiation in the j — ^ 1 
transition exhibits a nearly identical angular profile to the 
ji = 1 — >■ 2 transition. 

Figure [U] shows the level population of the excited state, 
ni/uo as a function of time in units of Tjef. Since the up- 
per level, j = 2, remains negligibly populated throughout 
the calculation, the ground state population is no /no ~ 
1 — ni/rio- The top panel shows ni/rio as a function of 
t/Trcf at several reference depths, r = 0,2,5. The bottom 
panel shows n\/no as a function of r for several times, 
t = Tic, l.STic, 2Tic. The population of the excited state 
peaks at t/Tiei ~ 500, which is when the radiation field close 
to the T — boundary has its maximum. The j = — > 2 
transition, followed by a spontaneous decay j = 2 — > 1, is the 
primary mechanism for populating the j = 1 level. Since the 
radiation field in the j = ^ 2 transition rapidly decreases 
at line centre for r ^ 1, most of the excitation events occur 
at T < 1. This explains why the population of the excited 
state peaks with the external radiation in the top panel: the 
light travel time to the layers in which significant excitation 
occurs is relatively small compared to the timescales over 



Figure 9. Level population of the first excited state, ni/rio as a 
function of time in units of T^ef , for several reference depths (top 
panel) ; and as a function of reference depth for several times (bot- 
tom panel) . All calculations have been performed for the canonical 
Model III of ^5:21 



which variation occurs. It is interesting to note that, as seen 
in the t = Tlc curve in the bottom panel of Figure [9] the 
population ni actually has its maximum value at r « 0.2. 
This is due to the fact that the external radiation has not 
yet been significantly extinguished, while the diffuse field is 
generated by absorption and re-emission into non-zero an- 
gles. The combination of these effects causes the radiation 
field in the j = — ^ 2 transition to peak at r > 0, leading 
to a maximum in the ni population at the same point. 

Figure[lO]shows the time development of the j = 1 — > 2, 
u„ line profile as a function of x, at r = Tmax with fi — 
1. The intensity scale is normalized to the maximum value 
of the radiation field, which occurs at t/Tid ~ 1500 for 
j: = 4. At line center (x = 0), the field decreases to its 
minimum intensity approximately when photons from the 
peak emission of the external source cross the medium; as 
expected, this is when the field in the line wing {x — 4) 
achieves its maximum intensity. When the medium contains 
a significant fraction of atoms in the j ~ 1 state, absorption 
at line center decreases the intensity of the radiation field at 
earlier times than for frequencies a; > 0. At t/T^cf > 1800, 
when the atoms have depopulated to the ground state, the 
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Figure 10. Time development of the line profile as a function 
of X. Results are shown for r = Tmax and fi = 1. The time is 
given in units of t/Tref, while the intensity scale is normalized 
to the maximum value of the radiation field, which occurs at 
i/Tvef ~ 1500 for x = 4. The calculations were performed for the 
canonical Model III of ^5^21 
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Figure 11. Line decrement at ^ = 1, r = Tmax, defined as the 
radiation variable Mj, at line centre divided by its value in the 
line wing {x = 4), Mccnt /"wing • The decrement is shown as a 
function of time in units of T^ef- The solid and dashed curves 
show results from the full calculation and coUimated approxima- 
tion, respectively. The dotted curve shows the relative difference 
between the full calculation and approximation. All calculations 
have been performed for the canonical Model III of i|5.2l 



time variation of photons at all frequencies is once again the 
same. 

Figure [11] shows ticcnt / Uwing for the j = 1 — >■ 2 transi- 
tion, at r = Tmax and /i = 1, as a function of time (solid 
curve). The decrement can be used to define an efi'ective 
optical depth, 



Tcfl 



' log(Uccnt/?iwing); 



(68) 



for the canonical Model III, t^s ~ 0.3. As expected, the peak 
effective optical depth in the j = 1 — ^ 2 transition occurs 
at time t/Trd ~ Tmax + 500, which corresponds to the time 



at which photons from the peak of the external source cross 
the medium. These photons experience large extinction at 
line centre due to excitation of the atoms into the j — 1 
level. 

Several previous works in the literature have used sim- 
plified treatments of radiative transfer to infer the state 
of atomic populations in absorbing material around GRBs 
from the effective optical depth in absorption lines of af- 
terglow specra (see ijS}. These works neglected the dif- 
fuse field in their solution to the radiative transfer equa- 
tion; we refer to this solution as the ' c oUimated approxi- 
mation' (see, e.g.. IWeeswiik et al.ll2007l : iD'Elia et al]l2009l : 
iRobinson et al.ll2010^ . In this case, the external radiation 
is propagated through the slab, but angularly redistributed 
photons are neglected in determining the values of the level 
populations. In this approximation, the calculation becomes 
much simpler, and the convergence of the level populations is 
almost immediate. The result is shown by the dashed curve 
in Figure 1111 Clearly, significant errors are introduced into 
the calculation of the decrement when the diffuse component 
of the radiation field is neglected. This is due to the fact that 
neglecting the diffuse field implies a lower excitation rate of 
atoms from the ground state to j = 1. If Tvr is of order 
Tref , then the line centre extinction for photons propagating 
through the medium will be increased when the diffuse field 
is taken into account, resulting in a larger decrement. 



5.3 Variation of Model Parameters (Models 
IV-VI) 

In this section, we explore the effects of changing the atomic 
physics; the transient amplitude. A; and the value of the 
background radiation, Jb. 

We performed a similar calculation to that described 
in Model III, except that we multiplied the Einstein coeffi- 
cients for the j = 1 ^ 2 transition by a factor of ten (Model 
IV). Figure [T5] is identical to Figure [S] except that it shows 
results for Model IV. In addition, the results in the top panel 
were computed for t = 0, 10, 100. The figure indicates that 
the effective optical depth for the j = 1 ^ 2 transition is 
much larger for Model IV than for Model III. An interesting 
feature in the top panel of Figure [12] is the slight decrease in 
the intensity around t/T^^i = t + 500, where the radiation 
field peaks in Figure [6] The new feature occurs because of 
the increase in the population of the j = 1 level, resulting in 
larger extinction. Most photons with frequencies near line 
centre of the j = 1 ^ 2 transition that pass through spatial 
point T f» 500 at time t/Tici ~ t-I-500 originated in the peak 
of the external radiation. However, the intensity exhibits an 
overall decrease due to the enhanced extinction. Such a fea- 
ture was absent from Figure |S] because the total extinction 
was too small to compensate for the increase in intensity 
from the external boundary. The same phenomenon is re- 
sponsible for the dip in intensity around t ~ 500 in the 
solid curve {t = Tlc) of the bottom panel. 

Figure [13] shows the angular distribution of radiation 
for Model IV. The results are similar to those in Figure [S] 
for Model III, with the radiation field exhibiting increased 
isotropy at t = Tmax, due to increased angular redistribution 
reflected in the larger value of the effective optical depth. 

Figure [14] is identical to the top panel of Figure [9] ex- 
cept that it shows ni/jio as a function of time at t = 0, 1, 5. 
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Figure 12. Same as Figure [B] except tiiat results are shown for 
Model IV of ESI 



Figure 13. Same as Figure [S] except that results are shown for 
Model IV of ^531 



At small reference depths, we see an increase in the maxi- 
mum population of level j = 1; this is due to the increase 
in the spontaneous emission coefficient for the j = 2 — >■ 1 
transition. 

Figure [T^ shows the decrement in the j = 1 — >■ 2 transi- 
tion as a function of time in units of Trcf. The plot exhibits 
the same overall features as Figure [TT] except that the maxi- 
mum effective optical depth at t/T^ci ~ Tniax+500 is T^g ~ 3. 
The larger Einstein coefficients for the j — 1 ^ 2 transition 
cause increased extinction per excitation event. Thus, while 
both the full calculation and the coUimated approximation 
exhibit enhanced extinction, there is a larger relative error 
associated with employing the latter (up to 100%) compared 
to Model III. This is due to the fact that there are fewer ex- 
citations to the j — 1 level in the coUimated approximation 
due to the neglect of the diffuse field, and this error is am- 
plified in the extinction by the change in atomic physics for 
Model IV. 

We performed another calculation using the same pa- 
rameters and Einstein coefficients as in Model III, except for 
an increase in the amplitude of the transient component of 
equation ([6| to A = 1 (Model V). 

Figure [in] is identical to Figure [T51 except that it shows 
results for Model V. The maximum effective optical depth 
in the j — 1 ^ 2 transition is TcH ~ 3. However, the rela- 




Figure 14. Same as the top panel of Figure [T4l except that results 
arc shown for Model IV of i|5.3l 
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Figure 15. Same as Figure [TTI except that results are shown for 
Model IV of 3531 



Figure 17. Same as Figure [TTI except that results are shown for 
Model VI of [ES] 
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Figure 16. Same as Figure [TTI except that results are shown for 
Model V of 



tive error introduced into the calculation by neglecting the 
diffuse field has increased by a factor of six. This is due to 
the increase in the maximum amplitude of the external ra- 
diation. In Model IV, the increased extinction relative to 
Model III resulted from the stronger j = 1 — >■ 2 transition 
in the atomic model, which affected both the full calcula- 
tion and the coUimated approximation. In Model V, the in- 
creased extinction relative to Model III is due to enhanced 
excitation to the j = 1 level; this affects the full calculation 
and the coUimated approximation differently. In the former, 
inclusion of the diffuse field fully incorporates the effect of 
increasing the external radiation amplitude. In the latter, 
a significant fraction of the amplitude increase is neglected, 
leading to a smaller effective optical depth for the coUimated 
approximation relative to Model IV. 

We explore a final case in which the background ra- 
diation is a significant fraction of the transient amplitude 
(Model VI). For this case, we change the parameters of the 
canonical Model III to Ib/Io = 0.1 and A/Io = 1. 

For the most part, the results for Model VI are similar to 



those for Model V. This is not unexpected as the maximum 
amplitude of the transient component of the external source 
is the same in both cases. One major difference is that the 
system begins and relaxes to an equilibrium state in which 
the diffuse field and excited j — 1 level population are non- 
zero. At reference depth r = for times t/Tici ~ and 
t/Tj-ci > 10^, the excited state populations are ni/rio ~ 0.08 
(note that this is approximately equal to the maximum pop- 
ulation at r = for Model III). Figure [TTI shows the line 
decrement results for Model VI. At t < Tlc, the line decre- 
ment is equal to its equilibrium value, itcont /wwing ~ 0.75. 
After a light crossing time, it decreases until the effective 
optical depth is Tch ~ 3, after which it relaxes back to its 
equilibrium value. 

The coUimated approximation was altered slightly for 
the calculations of Model VI. In this case, the steady-state 
diffuse field was determined at t/T^^i = and held fixed 
during the calculation; only the coUimated radiation was 
allowed to vary dynamically. As in Model V, the relative 
error introduced by neglecting changes in the diffuse field is 
significant, greater than 600% at maximum. 



6 DISCUSSION 

We have developed a general formalism for the numerical so- 
lution of the TDRTE and atomic rate equations using ALL 
We applied our method to a three-level atomic model that 
included allowed transitions between the ground and first ex- 
cited states to a common upper level, as well as a forbidden 
transition between the two former states. Our calculations 
exhibit convergence of the iterative ALI method and, when 
augmented with additional acceleration techniques, can be 
used in efficient calculations of line transfer with realistic 
atomic models. Our results show time coupled effects in the 
evolution of the atomic level populations, radiation field, 
and emergent spectra. We have compared our solution to 
an approximate treatment in which photons that are redis- 
tributed in angle are neglected (i.e., the diffuse field is ig- 
nored) , and showed that this approximation results in errors 
of 20 — 600% in the effective optical depth for the j = 1 2 
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transition. This error increases with the magnitude of Tcff . 
The effect occurs because the diffuse field photons affect the 
level populations and extinction coefficient for the excited 
state transition. If the time variability of the source falls 
in the intermediate regime (see then photons passing 
through the medium along the line of sight experience in- 
creased extinction when the diffuse radiation field is included 
in the calculation. 

These conclusions are relevant for timescales obeying 
the relations Tvr > Trcf and Tvr < Tlc- Assuming a 
Doppler line profile in complete redistribution, we can write 
the ratio Tva/Trcf as 



Trcf 



0.893 



rvR 

103 S 
Bo2 



fO~ 



b/c 



10^ cm^ s~i erg" 



-i) (cm-a) 



(69) 



and the maximum optical depth as 



1014cm 

Bq2 



10^ cm^ s~i erg 



6/c 



(70) 



If the timescales characterizing the external radiation source 
and medium satisfy the intermediate regime inequalities, 
then radiation transfer in the astrophysical system should 
be treated using the ALI methods developed in this paper. 
If, however, Tvr ^ Trcf, then the collimated approximation 
can be used without loss of accuracy. This is due to the fact 
that the diffuse radiation field is generated on timescales 
> Trcf, and if the transient pulse propagates through the 
medium on smaller timescales, then the level populations 
will not be altered by the diffuse field in time to affect ex- 
tinction of the external photons. Conversely, if Tvr ^ TLc, 
then the radiation field is determined by the steady-state 
RTE. 

Using the scaling of the quantities in equations (|69p and 
H70p . we can draw conclusions about what kind of astrophys- 
ical systems require the full TDRTE solution presented in 
this paper. In all cases, we assume that the characteristic 
transition has an Einstein coefficient _Bo2 ^ 10^ cm^ s~i 
erg"^ 

In absorbing systems composed of metal ions with Uo ~ 
10~^ cm"'^, if the width of the medium is W > f pc (corre- 
sponding to a minimum Tmax ~ 1), then the time variability 
of the external source must be Tvr > 10* s for the full ALI 
solution to apply. In contrast, for absorbing systems at a 
higher density, no ~ 1 cm~^, the full ALI solution applies 
when Tvr ~ 10^ s and > 3 x 10^^ cm. 

Time variability is common in astrophysical sources. 
Notable examples of highly transient emitters are GRBs, 
which emit prompt 7-ray emission at early times, followed by 
longer wavelength afterglow radiation in the X-ray, UV and 
optical bands. The afterglow emission occurs on timescales 
that range from minutes to days. A large number of stud- 
ies have examined the time dependence of afterglow ab- 
sorption spectra due to photoionization of the medium; 
this variability occurs on timescales c omparable to the ob- 
servation time ( Perna fc Loebl [l99i: 'Bottcher. et al ] I1999I: 
Perna fc Lazzati|[20 02: Laz zati fc Pern a 2002; Mirabal et aP 
2002l : lRobin'son' et al.ll2010l ). Observations of time-dependent 



absorption h a ve be e n reported for a number of GRBs 
jThone et al.l I2OIII : IVreeswiik et al.l |2012| : iDe Cia et all 
I2OI2I I. In a few cases, excitation of fine line transitions of 
Fe II we re observed in the abs e nce of appreciable photoion- 
ization (|Vreeswiik et all l2007l : iD'Elia et al.l |2009| ). So far, 
none of the time-dependent radiative transfer calculations 
developed specifically for GRBs has included source terms 
for the diffuse field. Previous work emphasized the modelling 
of time-dependent absorption effects due to photoionization, 
rather than the details of line transfer. 

Consider the results of our analysis for a GRB vari- 
ability timescale of TLc ~ lO'^ s, and typical interstellar 
medium hydrogen density of wh ~ 1 cm~"^. For metal ions 
with Uo ~ 10~^ cm~^, the condition Tvr ^ Tref is satisfied, 
and the collimated approximation is appropriate. However, 
for hydrogen, Tvr ~ Trcf, and the diffuse field can affect line 
transfer through the absorber, coupling to the photoioniza- 
tion state. Thus, more detailed studies, using the methods 
described in our paper, should be performed when interpret- 
ing observations of afterglow line variability. In addition, if 
the GRB occurs in a molecular cloud, with densities as high 
as ~ 10^ cm^'^, then a full solution to the time-dependent 
line transfer problem may be needed even for metal lines. 

AGNs form another important class of time variable 
sources. Typical AGN soft X-ray variability timescales range 
from ~ 1 ksec in nearby, small black hole mass S eyfert galajc- 
ies, t o ~ 100 ksec in distant, massive quasars |Fiore et al.l 
Il998l ). The central source is surrounded by a cloud of ab- 
sorbers with hydrogen densities as high as ~ 10^ — 10* cm. 



located at distances of ~ 10^'^ — 10^* cm (|Baldwin et al.l 
ll995l : lKrongold et al.|[2007h . If the source variability is about 
1 ksec for an absorber with density rio ^ 10'' cm~^ (in 
metal lines), then Tvr ^ lO^Tof. In this case, the full 
time-dependent formulation of line transfer is needed when 
Tmax 10^ , implying an absorber width of VK > 3 x 10^'' cm. 
If the source variability timescale is longer, with Tvr ~ 
100 ksec, then the full solution is required when Tmax ^ 10^, 
implying an absorber width of IV > 3 x 10^^ cm. 

The examples above illustrate that, for any time vari- 
able source, there is a range of absorber properties for which 
line transfer must be treated by the ALI method. 

In Appendix [X] we further generalize our formalism to 
include a broad class of Runge-Kutta methods for the time 
integrations employed in ^ These techniques will enable fu- 
ture calculations that monitor the error in the time grid dis- 
cretization, or implement an adaptive time step algorithm. 
Such improvements can improve the stability of the radi- 
ation field integration, reduce artificial oscillations in the 
solution, and allow more efficient calculations for realistic 
atomic models. 
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APPENDIX A: RUNGE-KUTTA METHODS 

In this Appendix, we outline how the iterative scheme of 
tj4.4l can be impleme nted with implicit Runge-Kutta meth- 
ods (see, e.g., §11.7 of iHairer et al.lll99ll ). For simplicity, we 
neglect collisions in the rate equations and background con- 
tributions to the emissivity. Including these terms requires 
a straightforward generalization of the work below. 

Applying the method of lines to the simultaneous sys- 
tem of equations yields (c.f., i|4.2|) : 

.fj,d = '^{Rij.d.ni^d - Rji,dnj,d), (Al) 
I 

fu,di = CXdi («d+l/2,i - t'd-l/2,i) /Ardj 

L (A2) 

— Udi + Sdi , 
fv,d+l/2,i = CXd+l/2,i (Md+l,i ~Ud,) / Ar^ 

L (A3) 

— i^ti+l/2,i ■ 

Note that we used a second order spatial discretization 
above, though more complicated methods are possible. 
Higher order methods will lead to a more complicated spa- 
tial structure for the linear systems of equations (see below) . 

The Runge-Kutta solution for equation ()Aip can be 
written 

nj^kd = nj^k-i,d 

p-i 

+ Atk bpfj^d {tk-i + CpAtk,nl^^^) , (A4) 

with 

p-i 

^^,kd ~ ^j,k-l,d + Atk a-pp' fj,d(^tk-l 

p'=0 

+ Cp^*fe''^j,M'<,'fcd.---). (A5) 

where P is the number of stages, and Cp, 
the coefficients for the Runge-Kutta method. The nota- 
tion kj,nf . . . indicates that the p'**^ stage quantities 
for each level are substituted in the expression for fj^d- In 
equations (|A4|) and (|A5|) . contains a dependence on Ud 
through the quantities Rid,d defined in 

We can use the Runge-Kutta solution to equation HA3|I 
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to simplify the dependence on v in (|A2|) : 

p-i 

'"k,d+l/2,i = '"fc-l,d+l/2,i + Atfc ^ "pp"=Xfc,d+l/2,i 

p'=0 

X - - <d+i/2 J - (A6) 

where the quantities Arj^^f and ti+i/2 i evaluated at 

+ CpiAtk using the values This is a linear system 

of dimension P in the quantities ^1+1/2 r K w6 define the 
column vector 

_ / p=o p=p-i f\'7\ 

Vfc,d+l/2,i - (^«fc,d+l/2,i. • • •'"fc,d+l/2,ij ' 

we can write the system as 

Lfc,d+l/2,i ■ Vfc,d+l/2,i = 

where the elements of the P x P matrices hk,d+i/2.i and 
^k,d+i/2,i are determined by equation (|A6|) . The elements 
of the vector Vfc_i_d+i/2,i are all equal to Ufc_i,ti+i/2,i- The 
solution to equation (|A8[) can be obtained in 0{P^) opera- 
tions for each triple of indices (fc, d + 1/2, i). We write this 
solution as: 

Vfc,d+l/2,i ~ Aifc,d+l/2,i 

■ (Ufc,d+i_i - Ukdi) + A/'fe,d+l/2,i- (A9) 

We can substitute equation (jA9|) into the Runga-Kutta 
formula for u to obtain a formal solution for the diffuse radi- 
ation field. The Runge-Kutta solution for u can be written: 

p-i 

= Uk-i,di + Atk flpp'CxLi 

X [(*^fc,d-|-l/2,i ~ '"fc,d-l/2,i) / ^'''kdi 

-<d. + st]^ (AlO) 

Lfcdi • Ufedi = Mfcdi • (Vfc_£j+l/2,i — Vfc_d_i/2,i) 

-f Nfcdi • Sfedi + Ufc_i,di, (All) 

where the matrices and vectors are defined analogously to 
those in the equation for the v variable. If we substitute 
equation (|A9I) into HA11|) . we obtain a system of equations 
of the form: 

— Afcdi • Ufc,d_l,i -I- (Lkdi + A-kdi + Ckdi) ' Ukdi (A12) 
— Cfedi • Ufc_d+l,i ~ R-fcdi • Skdi + Bfcdi, 

where the vector Jikdi contains the dependences on the radi- 
ation variables at the previous time grid point, k — 1. Equa- 
tion (|A12|) can be supplemented by boundary conditions 
connecting the points d = 0, 1 and d = D — 2, D — 1. These 
conditions can be derived as in ^4.21 the method of lines is 
applied to first order, and the time integration proceeds as 
described above. Thus, we obtain a block tridiagonal linear 
system that can be solved in O (P'^D) operations for each 
(fc, i). We write the solution to equation (|A12|) as 

Uki — Afei ■ Ski + Bfci, (A13) 

where Aki is a D x D block matrix with P x P matrices as 



elements, and Jiki is a D x 1 block vector with P x 1 vectors 
as elements. 

To obtain a formal solution for the diffuse field, we aug- 
ment our time grid with the points {tk-i + CpAtk}- The 
solution then consists of the values u^^i ' ■ ■ ■ '^^d^^^ ' ^^^r 
k = 1, . . . K — I, where the quantities u^di^^^ represent the 
solution on the original K — 1 time grid points. Using the 
initial conditions, we can solve equation HA13|) successively 
for k — 1, . . . K — 1. The structure of this formal solution 
in terms of the source functions S^^ is quite complicated, 
coupling the various Runge-Kutta stages and spatial grid 
points. However, if we consider the ALI form of the formal 
solution, 

^Li = Afedi ('S'Li ~ ^kdi) + ^fcdi' (A14) 

we achieve a significant simplification. In equation (|A14|) . 
^kdi denotes the diagonal elements of the P x P matrices 
that are in turn along the diagonal of the DxD block matrix 
Afc . The source functions are evaluated with the level popu- 
lations kd' kd' ■ • • ' t'^^ current and previous (f) iter- 
ations. The diagonal elements of Aki can be obtained simul- 
taneously with the formal solution, usi ng the block matrix 
form o f the equations in Appendix B of lRvbicki fc Hummed 
As in ti4.4l we use the alternate form of the formal 
solution, 

'^fedi ~ '^kdi (jlkdi ~ Vkdi^ + ^kdi 

= ffd. E Ui,- (<.*d - <kd) + <L - (A15) 
IJ 

where M^^, = AZJx^,- 

We can substitute equation ([AISP into equation HA5|I 
and precondition the resulting non-linear expressions as de- 
scribed in i]4.3l and i]4.4l The ALI method then proceeds as 
follows: starting with the initial conditions, Uk^o,di, nj^k=o,d, 
we solve for the formal solution at A; = 1, . . . , A" — 1. At each 
k, we iterate over the populations ^.^ for all the atomic 
levels and Runge-Kutta stages, using the formal solution 
for u^Ij, approximate matrix A^^, and ALI linear system. 
The method is completely analogous to the technique de- 
scribed in i^4.4l indeed the equations in the main text are 
special cases of Runge-Kutta integrations with coefficients 

1 1 

given by the Butcher tableaux (Backward Euler) 

1 


and 1 1/2 (Crank-Nicholson). 

1/2 1/2 



